Imports System.Drawing
Imports System.Collections.Generic
Namespace Numerical
  Public Class Geometry
    ''' <summary>
    ''' Compute the intersection between the line from pt1 to pt2 and pt3 to pt4
    ''' </summary>
    ''' <param name="pt1"></param>
    ''' <param name="pt2"></param>
    ''' <param name="pt3"></param>
    ''' <param name="pt4"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Shared Function PolygonIntersection(ByVal pt1 As PointF, ByVal pt2 As PointF, ByVal pt3 As PointF, ByVal pt4 As PointF, ByRef intersectionPt As PointF) As Boolean
      Dim deltaY1 As Single = pt2.Y - pt1.Y
      Dim deltaX1 As Single = pt2.X - pt1.X
      Dim deltaY2 As Single = pt4.Y - pt3.Y
      Dim deltaX2 As Single = pt4.X - pt3.X
      Dim denominator As Single = deltaY2 * deltaX1 - deltaX2 * deltaY1

      If denominator = 0 Then
        'No intersection, parallel lines
        Return False
      Else
        Dim numerator As Single = deltaX2 * (pt1.Y - pt3.Y) - deltaY2 * (pt1.X - pt3.X)
        Dim t As Single = numerator / denominator
        If t >= 0 And t <= 1 Then
          intersectionPt.X = pt1.X + t * deltaX1
          intersectionPt.Y = pt1.Y + t * deltaY1
          Return True
        Else
          Return False
        End If
      End If
    End Function

    ''' <summary>
    ''' Determine the centroid of a series of points, with the area already precomputed
    ''' </summary>
    ''' <param name="sArea">The already computed area</param>
    ''' <returns>A <b>pointf</b> with the centroid</returns>
    ''' <remarks>
    ''' </remarks>
    Public Shared Function PolygonCentroid(ByVal pts As IList(Of PointF), ByVal sArea As Single) As PointF
      Dim i As Integer, oCentroid As PointF

      For i = 0 To pts.Count - 2
        Dim iTemp As Single
        iTemp = pts(i).X * pts(i + 1).Y - pts(i + 1).X * pts(i).Y
        oCentroid.X += (pts(i).X + pts(i + 1).X) * iTemp
        oCentroid.Y += (pts(i).Y + pts(i + 1).Y) * iTemp
      Next

      If sArea > 0 Then
        oCentroid.X = CInt(Math.Round(Math.Abs(oCentroid.X / (6 * sArea))))
        oCentroid.Y = CInt(Math.Round(Math.Abs(oCentroid.Y / (6 * sArea))))
      Else
        oCentroid.X = 0
        oCentroid.Y = 0
      End If

      Debug.WriteLine("Centroid is " & oCentroid.ToString, "Geometry.PolygonCentroid")
      Return oCentroid
    End Function

    ''' <summary>
    ''' Determine the area taken up by the polygon, in pixels units
    ''' </summary>
    ''' <returns>The area</returns>
    ''' <remarks>
    ''' </remarks>
    Public Shared Function PolygonArea(ByVal pts As IList(Of PointF)) As Double
      Dim i As Integer, sArea As Single

      If pts.Count > 2 Then
        For i = 0 To pts.Count - 2
          sArea += pts(i).X * pts(i + 1).Y - pts(i + 1).X * pts(i).Y
        Next
      End If

      sArea = Math.Abs(sArea / 2)
      Debug.WriteLine("Area is " & sArea, "Geometry.PolygonArea")
      Return sArea
    End Function

    ''' <summary>
    ''' Determine the border length, in pixel units.
    ''' </summary>
    ''' <returns>The length</returns>
    ''' <remarks>
    ''' </remarks>
    Public Shared Function PolygonPerimeter(ByVal pts As IList(Of PointF)) As Double
      'Compute the length of annotation in pixel units
      Dim sLength As Double, i As Integer
      For i = 1 To pts.Count - 1
        sLength += Math.Sqrt(Math.Pow(pts(i).X - pts(i - 1).X, 2) + Math.Pow(pts(i).Y - pts(i - 1).Y, 2))
      Next
      Debug.WriteLine("Perimeter is " & sLength, "AnnotationPoints.GetLength")
      Return sLength
    End Function

    ''' <summary>
    ''' Return the shape compactness parameter, a value between 0 - 1. A circle has 1, a square 0.9. 
    ''' A rectangle with one dimension twice the other has compactness 0.83.
    ''' Use this version when the area and perimeter have already been computed
    ''' </summary>
    ''' <param name="perimeter"></param>
    ''' <param name="area"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Shared Function PolygonCompactness(ByVal pts As IList(Of PointF), ByVal perimeter As Double, ByVal area As Double) As Double
      Dim compact As Double = Math.Sqrt(area) * 3.53 / perimeter
      Debug.WriteLine("Polygon has compactness " & compact)
      Return (compact)
    End Function

    ''' <summary>
    ''' Return the shape compactness parameter, a value between 0 - 1. A circle has 1, a square 0.9. 
    ''' A rectangle with one dimension twice the other has compactness 0.83.
    ''' </summary>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Shared Function PolygonCompactness(ByVal pts As IList(Of PointF)) As Double
      Return PolygonCompactness(pts, PolygonPerimeter(pts), PolygonArea(pts))
    End Function

    ''' <summary>
    ''' Return the bisector of an angle defined by three points. Result is a normalized vector from the origin.
    ''' </summary>
    ''' <param name="pt1"></param>
    ''' <param name="pt2"></param>
    ''' <param name="pt3"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Shared Function Bisector(ByVal pt1 As PointF, ByVal pt2 As PointF, ByVal pt3 As PointF) As PointF
      Dim v1 As New PointF(pt3.X - pt1.X, pt3.Y - pt1.Y)
      Dim d1, d2 As Double
      d1 = Algebra.Norm(pt2.X - pt1.X, pt2.Y - pt1.Y)
      d2 = Algebra.Norm(pt3.X - pt2.X, pt3.Y - pt2.Y)
      If d1 > 0 And d2 > 0 Then
        Dim ratio As Double = d1 / (d1 + d2)
        Dim bs As New PointF(pt1.X + ratio * v1.X - pt2.X, pt1.Y + ratio * v1.Y - pt2.Y)
        Dim d3 As Double = Algebra.Norm(bs.X, bs.Y)
        If d3 > 0 Then
          bs.X /= d3
          bs.Y /= d3
        Else
          'Points are collinear.
          bs.X = (pt1.Y - pt2.Y) / d1
          bs.Y = (pt2.X - pt1.X) / d2
        End If
        Return (bs)
      Else
        Return New PointF(0, 0)
      End If
    End Function


    ''' <summary>
    ''' Determine if a point is inside the area enclosed by the Annotation. 
    ''' Instead of using a region like we do when testing for a point with integer coordinates, 
    ''' we use a point-in-polygon method (ref. W. F. Randolph). This can be much slower, especially when
    ''' there are many points in the path, but it will be much more precise and not suffer from 
    ''' rounding errors that are inherent when converting an annotation to a region.
    ''' </summary>
    ''' <param name="pt"></param>
    ''' <param name="pts"></param>
    ''' <returns></returns>
    ''' <remarks>IMPORTANT: Do not use the GraphicsPath.IsVisible method in .NET, 
    ''' it is buggy and does not work for floating point coordinates!</remarks>
    Public Shared Function PointInPolygon(ByVal pt As PointF, ByVal pts As IList(Of PointF)) As Boolean
      Dim pt1, pt2 As PointF
      Dim nrPts As Integer = pts.Count
      Dim result As Boolean = False
      For i As Integer = 0 To nrPts - 1
        pt1 = pts(i)
        If i = 0 Then
          pt2 = pts(nrPts - 1)
        Else
          pt2 = pts(i - 1)
        End If

        If ((((pt1.Y <= pt.Y) AndAlso (pt.Y < pt2.Y)) OrElse ((pt2.Y <= pt.Y) AndAlso (pt.Y < pt1.Y))) _
          AndAlso (pt.X < (pt2.X - pt1.X) * (pt.Y - pt1.Y) / (pt2.Y - pt1.Y) + pt1.X)) Then
          result = Not result
        End If
      Next
      Return result
    End Function

    ''' <summary>
    ''' Determine the location of a point versus a line defined by 2 points. If the result is 0 then the point is on the line. 
    ''' For result 1 it is in the positive halfspace to the right of the vector P1P2. For result -1 the point is in the 
    ''' negative halfspace to the left of the vector P1P2
    ''' </summary>
    ''' <param name="pt"></param>
    ''' <param name="pt1"></param>
    ''' <param name="pt2"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Shared Function HalfSpace(ByVal pt As PointF, ByVal pt1 As PointF, ByVal pt2 As PointF) As Double
      Return Math.Sign((pt2.Y - pt1.Y) * (pt.X - pt1.X) - (pt.Y - pt1.Y) * (pt2.X - pt1.X))
    End Function

  End Class
End Namespace


